########################################################################################################################################################
# April 22 2024
# Results Replication File
# Rebecca Cordell, Reed M. Wood and Thorin M. Wright
# Disease and Dissent: Epidemics as a Catalyst for Social Unrest
# Global Studies Quarterly
# https://academic.oup.com/isagsq
########################################################################################################################################################

# Clear work environment
rm(list=ls())

# Required Packages
require(glmmTMB)
require(texreg)
require(stats)
require(emmeans)
require(ggplot2)
require(cem)
options(scipen=999)

###############
# Read in data
###############
disease_dissent<-read.csv("CordellWoodWright_2023_DiseaseandDissent.csv", header= TRUE, stringsAsFactors = FALSE)

############################################################################
# Table 1. Multi-level Poisson regression results, epidemic binary variable
############################################################################

# Pre-matched sample
pre_matched<-disease_dissent

# Post-matched sample
post_matched<-disease_dissent[disease_dissent$matched==1,]

# Model 1 Pre-matched sample
tab1_m1 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = poisson, data=pre_matched)

# Model 2 Pre-matched sample
tab1_m2 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = poisson, data=pre_matched)

# Model 3 Post-matched sample
tab1_m3 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = poisson, data=post_matched)

# Model 4 Post-matched sample
tab1_m4 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = poisson, data=post_matched)

# Extract results
source(system.file("other_methods","extract.R",package="glmmTMB"))
extract.tab1_m1<-extract.glmmTMB(tab1_m1)
extract.tab1_m2<-extract.glmmTMB(tab1_m2)
extract.tab1_m3<-extract.glmmTMB(tab1_m3)
extract.tab1_m4<-extract.glmmTMB(tab1_m4)

# Print Regression Table
screenreg(list(extract.tab1_m1, extract.tab1_m2, extract.tab1_m3, extract.tab1_m4), custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"), stars = c(0.01, 0.05, 0.1))

# Observations
n_tab1_m1<-nrow(tab1_m1$frame)
n_tab1_m2<-nrow(tab1_m2$frame)
n_tab1_m3<-nrow(tab1_m3$frame)
n_tab1_m4<-nrow(tab1_m4$frame)

# Log Likelihood
log_tab1_m1<-logLik(tab1_m1)
log_tab1_m2<-logLik(tab1_m2)
log_tab1_m3<-logLik(tab1_m3)
log_tab1_m4<-logLik(tab1_m4)

# AIC
aic_tab1_m1<-AIC(tab1_m1)
aic_tab1_m2<-AIC(tab1_m2)
aic_tab1_m3<-AIC(tab1_m3)
aic_tab1_m4<-AIC(tab1_m4)

# BIC
bic_tab1_m1<-BIC(tab1_m1)
bic_tab1_m2<-BIC(tab1_m2)
bic_tab1_m3<-BIC(tab1_m3)
bic_tab1_m4<-BIC(tab1_m4)

# Create Model Fit Stats Table
all_tab1_m1<-rbind(n_tab1_m1,log_tab1_m1,aic_tab1_m1,bic_tab1_m1)
all_tab1_m2<-rbind(n_tab1_m2,log_tab1_m2,aic_tab1_m2,bic_tab1_m2)
all_tab1_m3<-rbind(n_tab1_m3,log_tab1_m3,aic_tab1_m3,bic_tab1_m3)
all_tab1_m4<-rbind(n_tab1_m4,log_tab1_m4,aic_tab1_m4,bic_tab1_m4)

# Print Model Fit Stats
all_tab1<-round(cbind(all_tab1_m1,all_tab1_m2,all_tab1_m3,all_tab1_m4),digits=2)
all_tab1

##############################################################################
# Table 2. Multi-level Poisson regression results, epidemic severity variable
##############################################################################

# Pre-matched sample
pre_matched$epidemic_severity<-as.factor(pre_matched$epidemic_severity)

# Post-matched sample
post_matched$epidemic_severity<-as.factor(post_matched$epidemic_severity)

# Model 1 Pre-matched sample
tab2_m1 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = poisson, data=pre_matched)

# Model 2 Pre-matched sample
tab2_m2 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = poisson, data=pre_matched)

# Model 3 Post-matched sample
tab2_m3 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = poisson, data=post_matched)

# Model 4 Post-matched sample
tab2_m4 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = poisson, data=post_matched)

# Extract results
source(system.file("other_methods","extract.R",package="glmmTMB"))
extract.tab2_m1<-extract.glmmTMB(tab2_m1)
extract.tab2_m2<-extract.glmmTMB(tab2_m2)
extract.tab2_m3<-extract.glmmTMB(tab2_m3)
extract.tab2_m4<-extract.glmmTMB(tab2_m4)

# Print Regression Table
screenreg(list(extract.tab2_m1, extract.tab2_m2, extract.tab2_m3, extract.tab2_m4), custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"), stars = c(0.01, 0.05, 0.1))

# Observations
n_tab2_m1<-nrow(tab2_m1$frame)
n_tab2_m2<-nrow(tab2_m2$frame)
n_tab2_m3<-nrow(tab2_m3$frame)
n_tab2_m4<-nrow(tab2_m4$frame)

# Log Likelihood
log_tab2_m1<-logLik(tab2_m1)
log_tab2_m2<-logLik(tab2_m2)
log_tab2_m3<-logLik(tab2_m3)
log_tab2_m4<-logLik(tab2_m4)

# AIC
aic_tab2_m1<-AIC(tab2_m1)
aic_tab2_m2<-AIC(tab2_m2)
aic_tab2_m3<-AIC(tab2_m3)
aic_tab2_m4<-AIC(tab2_m4)

# BIC
bic_tab2_m1<-BIC(tab2_m1)
bic_tab2_m2<-BIC(tab2_m2)
bic_tab2_m3<-BIC(tab2_m3)
bic_tab2_m4<-BIC(tab2_m4)

# Create Model Fit Stats Table
all_tab2_m1<-rbind(n_tab2_m1,log_tab2_m1,aic_tab2_m1,bic_tab2_m1)
all_tab2_m2<-rbind(n_tab2_m2,log_tab2_m2,aic_tab2_m2,bic_tab2_m2)
all_tab2_m3<-rbind(n_tab2_m3,log_tab2_m3,aic_tab2_m3,bic_tab2_m3)
all_tab2_m4<-rbind(n_tab2_m4,log_tab2_m4,aic_tab2_m4,bic_tab2_m4)

# Print Model Fit Stats
all_tab2<-round(cbind(all_tab2_m1,all_tab2_m2,all_tab2_m3,all_tab2_m4),digits=2)
all_tab2

#############################################################################################################################################################################################
# P.8 "For our Epidemic binary variable, the probability of observing a social unrest event increases by 18 percent when a first-order administrative unit experiences an epidemic outbreak." 
#############################################################################################################################################################################################

# Calculate Marginal Effects
dummy.me<-emmeans(tab1_m4, specs = pairwise ~ epidemic_dummy)
dummy.me<- eff_size(dummy.me, sigma(tab1_m4), edf = df.residual(tab1_m4))
dummy.me<-as.data.frame(dummy.me)
dummy.me$effect.size <- gsub("-", "", dummy.me$effect.size)
dummy.me$effect.size<-as.numeric(dummy.me$effect.size)

# Print
round(dummy.me$effect.size, digits=2)*100

#########################################################################################################################################################################################################################################################################################################################
# P.8 "Based on the marginal effects presented in figure 3, the probability of observing a single Social Unrest event is expected to increase by 17 percent when a first-order administrative unit experiences an epidemic outbreak of low severity and 18 percent when it experiences an epidemic of moderate severity." 
#########################################################################################################################################################################################################################################################################################################################

# Calculate Marginal Effects
severity.me<-emmeans(tab2_m4, specs = pairwise ~ epidemic_severity)
severity.me<- eff_size(severity.me, sigma(tab2_m4), edf = df.residual(tab2_m4))
severity.me<-as.data.frame(severity.me)

# Create Marginal Effects Table for Figure 3
severity.me_1_2<-severity.me[1:2,]
severity.me_1_2$effect.size <- gsub("-", "", severity.me_1_2$effect.size)
severity.me_1_2$lower.CL <- gsub("-*", "", severity.me_1_2$lower.CL)
severity.me_1_2$upper.CL <- gsub("-", "", severity.me_1_2$upper.CL)
severity.me_3<-severity.me[3,]
severity.me_3 <- data.frame("(0 - 3)","-0.119842027118805", "0.11257795", "3722","0.100878485763471", "-0.340562540001081")      
names(severity.me_3) <- c("contrast", "effect.size", "SE", "df", "lower.CL", "upper.CL") 
severity.me<-rbind(severity.me_1_2, severity.me_3)
severity.me$effect.size<-as.numeric(severity.me$effect.size)
severity.me$lower.CL<-as.numeric(severity.me$lower.CL)
severity.me$upper.CL<-as.numeric(severity.me$upper.CL)
severity.me$effect.size<-as.numeric(severity.me$effect.size)

# Print 
round(severity.me$effect.size[1], digits=2)*100

##############################################################################
# Figure 3. Average marginal effect of epidemic severity variable
##############################################################################

# Print Figure 3
ggplot(severity.me, aes(contrast, effect.size)) +
  geom_point() +
  geom_errorbar(aes(ymin=lower.CL, ymax=upper.CL), width=.1,
                position=position_dodge(.9)) +
  xlab("") +
  ylab("Average Marginal Effects") +
  scale_x_discrete(labels=c("(0 - 1)" = "Low Severity",
                            "(0 - 2)" = "Moderate Severity", 
                            "(0 - 3)" = "Severe")) +
  theme_bw() +
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL), alpha = .1) +
  geom_hline(yintercept = 0, linewidth=0.1)

###########################################################################
# Table A1
# Descriptive Statistics of Independent and Control Variables
###########################################################################

data<-disease_dissent
data$epidemic_severity<-as.character(data$epidemic_severity)
data$epidemic_severity<-as.numeric(data$epidemic_severity)

# N
epidemic_dummy <- data[!is.na(data$epidemic_dummy),]
epidemic_severity <- data[!is.na(data$epidemic_severity),]
pop_dens_log <- data[!is.na(data$pop_dens_log),]
night_lights_log <- data[!is.na(data$night_lights_log),]
conflict_intens_log <- data[!is.na(data$conflict_intens_log),]
democracy <- data[!is.na(data$democracy),]
gdppc_log <- data[!is.na(data$gdppc_log),]
pop_log <- data[!is.na(data$pop_log),]
federalism <- data[!is.na(data$federalism),]
hiv_log <- data[!is.na(data$hiv_log),]
temporal_lag_social_unrest <- data[!is.na(data$temporal_lag_social_unrest),]
temporal_lag_repression  <- data[!is.na(data$temporal_lag_repression ),]
spatial_lag_social_unrest  <- data[!is.na(data$spatial_lag_social_unrest ),]
political_exclud_group  <- data[!is.na(data$political_exclud_group ),]
gdp_growth  <- data[!is.na(data$gdp_growth ),]
press_freedom <-data[!is.na(data$press_freedom),]

n_epidemic_dummy<-length(epidemic_dummy$epidemic_dummy)
n_epidemic_severity<-length(epidemic_severity$epidemic_severity)
n_pop_dens_log<-length(pop_dens_log$pop_dens_log)
n_night_lights_log<-length(night_lights_log$night_lights_log)
n_conflict_intens_log<-length(conflict_intens_log$conflict_intens_log)
n_democracy<-length(democracy$democracy)
n_gdppc_log<-length(gdppc_log$gdppc_log)
n_pop_log<-length(pop_log$pop_log)
n_federalism<-length(federalism$federalism)
n_hiv_log<-length(hiv_log$hiv_log)
n_temporal_lag_social_unrest<-length(temporal_lag_social_unrest$temporal_lag_social_unrest)
n_temporal_lag_repression <-length(temporal_lag_repression $temporal_lag_repression )
n_spatial_lag_social_unrest <-length(spatial_lag_social_unrest $spatial_lag_social_unrest )
n_political_exclud_group <-length(political_exclud_group$political_exclud_group  )
n_gdp_growth <-length(gdp_growth$gdp_growth)
n_press_freedom <-length(press_freedom$press_freedom)

# Standard Deviation
sd_epidemic_dummy<-sd(epidemic_dummy$epidemic_dummy) 
sd_epidemic_severity<-sd(epidemic_severity$epidemic_severity) 
sd_pop_dens_log<-sd(pop_dens_log$pop_dens_log)
sd_night_lights_log<-sd(night_lights_log$night_lights_log)
sd_conflict_intens_log<-sd(conflict_intens_log$conflict_intens_log)
sd_democracy<-sd(democracy$democracy)
sd_gdppc_log<-sd(gdppc_log$gdppc_log)
sd_pop_log<-sd(pop_log$pop_log)
sd_federalism<-sd(federalism$federalism)
sd_hiv_log<-sd(hiv_log$hiv_log) 
sd_temporal_lag_social_unrest<-sd(temporal_lag_social_unrest$temporal_lag_social_unrest)
sd_temporal_lag_repression <-sd(temporal_lag_repression $temporal_lag_repression )
sd_spatial_lag_social_unrest <-sd(spatial_lag_social_unrest $spatial_lag_social_unrest )
sd_political_exclud_group <-sd(political_exclud_group$political_exclud_group)
sd_gdp_growth <-sd(gdp_growth$gdp_growth)
sd_press_freedom <-sd(press_freedom$press_freedom)

# Mean
mean_epidemic_dummy<-mean(epidemic_dummy$epidemic_dummy) 
mean_epidemic_severity<-mean(epidemic_severity$epidemic_severity) 
mean_pop_dens_log<-mean(pop_dens_log$pop_dens_log)
mean_night_lights_log<-mean(night_lights_log$night_lights_log)
mean_conflict_intens_log<-mean(conflict_intens_log$conflict_intens_log)
mean_democracy<-mean(democracy$democracy)
mean_gdppc_log<-mean(gdppc_log$gdppc_log)
mean_pop_log<-mean(pop_log$pop_log)
mean_federalism<-mean(federalism$federalism)
mean_hiv_log<-mean(hiv_log$hiv_log) 
mean_temporal_lag_social_unrest<-mean(temporal_lag_social_unrest$temporal_lag_social_unrest) 
mean_temporal_lag_repression <-mean(temporal_lag_repression $temporal_lag_repression ) 
mean_spatial_lag_social_unrest <-mean(spatial_lag_social_unrest $spatial_lag_social_unrest ) 
mean_political_exclud_group <-mean(political_exclud_group$political_exclud_group) 
mean_gdp_growth <-mean(gdp_growth $gdp_growth) 
mean_press_freedom <-mean(press_freedom $press_freedom) 

# Min
min_epidemic_dummy<-min(epidemic_dummy$epidemic_dummy) 
min_epidemic_severity<-min(epidemic_severity$epidemic_severity) 
min_pop_dens_log<-min(pop_dens_log$pop_dens_log)
min_night_lights_log<-min(night_lights_log$night_lights_log)
min_conflict_intens_log<-min(conflict_intens_log$conflict_intens_log)
min_democracy<-min(democracy$democracy)
min_gdppc_log<-min(gdppc_log$gdppc_log)
min_pop_log<-min(pop_log$pop_log)
min_federalism<-min(federalism$federalism)
min_hiv_log<-min(hiv_log$hiv_log) 
min_temporal_lag_social_unrest<-min(temporal_lag_social_unrest$temporal_lag_social_unrest) 
min_temporal_lag_repression <-min(temporal_lag_repression $temporal_lag_repression ) 
min_spatial_lag_social_unrest <-min(spatial_lag_social_unrest $spatial_lag_social_unrest ) 
min_political_exclud_group <-min(political_exclud_group$political_exclud_group) 
min_gdp_growth <-min(gdp_growth$gdp_growth) 
min_press_freedom <-min(press_freedom$press_freedom) 

# Max
max_epidemic_dummy<-max(epidemic_dummy$epidemic_dummy) 
max_epidemic_severity<-max(epidemic_severity$epidemic_severity) 
max_pop_dens_log<-max(pop_dens_log$pop_dens_log)
max_night_lights_log<-max(night_lights_log$night_lights_log)
max_conflict_intens_log<-max(conflict_intens_log$conflict_intens_log)
max_democracy<-max(democracy$democracy)
max_gdppc_log<-max(gdppc_log$gdppc_log)
max_pop_log<-max(pop_log$pop_log)
max_federalism<-max(federalism$federalism)
max_hiv_log<-max(hiv_log$hiv_log) 
max_temporal_lag_social_unrest<-max(temporal_lag_social_unrest$temporal_lag_social_unrest) 
max_temporal_lag_repression <-max(temporal_lag_repression $temporal_lag_repression ) 
max_spatial_lag_social_unrest <-max(spatial_lag_social_unrest $spatial_lag_social_unrest ) 
max_political_exclud_group <-max(political_exclud_group$political_exclud_group) 
max_gdp_growth <-max(gdp_growth $gdp_growth) 
max_press_freedom <-max(press_freedom $press_freedom) 

# Create table
sumstats_epidemic_dummy<-cbind(n_epidemic_dummy, mean_epidemic_dummy, sd_epidemic_dummy, min_epidemic_dummy, max_epidemic_dummy)
sumstats_epidemic_severity<-cbind(n_epidemic_severity, mean_epidemic_severity, sd_epidemic_severity, min_epidemic_severity, max_epidemic_severity)
sumstats_pop_dens_log<-cbind(n_pop_dens_log, mean_pop_dens_log, sd_pop_dens_log, min_pop_dens_log, max_pop_dens_log)
sumstats_night_lights_log<-cbind(n_night_lights_log, mean_night_lights_log, sd_night_lights_log, min_night_lights_log, max_night_lights_log)
sumstats_conflict_intens_log<-cbind(n_conflict_intens_log, mean_conflict_intens_log, sd_conflict_intens_log, min_conflict_intens_log, max_conflict_intens_log)
sumstats_democracy<-cbind(n_democracy, mean_democracy, sd_democracy, min_democracy, max_democracy)
sumstats_gdppc_log<-cbind(n_gdppc_log, mean_gdppc_log, sd_gdppc_log, min_gdppc_log, max_gdppc_log)
sumstats_pop_log<-cbind(n_pop_log, mean_pop_log, sd_pop_log, min_pop_log, max_pop_log)
sumstats_federalism<-cbind(n_federalism, mean_federalism, sd_federalism, min_federalism, max_federalism)
sumstats_hiv_log<-cbind(n_hiv_log, mean_hiv_log, sd_hiv_log, min_hiv_log, max_hiv_log)
sumstats_temporal_lag_social_unrest<-cbind(n_temporal_lag_social_unrest, mean_temporal_lag_social_unrest, sd_temporal_lag_social_unrest, min_temporal_lag_social_unrest, max_temporal_lag_social_unrest)
sumstats_temporal_lag_repression <-cbind(n_temporal_lag_repression , mean_temporal_lag_repression , sd_temporal_lag_repression , min_temporal_lag_repression , max_temporal_lag_repression )
sumstats_spatial_lag_social_unrest <-cbind(n_spatial_lag_social_unrest , mean_spatial_lag_social_unrest , sd_spatial_lag_social_unrest , min_spatial_lag_social_unrest , max_spatial_lag_social_unrest )
sumstats_political_exclud_group <-cbind(n_political_exclud_group , mean_political_exclud_group , sd_political_exclud_group , min_political_exclud_group , max_political_exclud_group )
sumstats_gdp_growth <-cbind(n_gdp_growth , mean_gdp_growth , sd_gdp_growth , min_gdp_growth , max_gdp_growth )
sumstats_press_freedom <-cbind(n_press_freedom , mean_press_freedom , sd_press_freedom , min_press_freedom , max_press_freedom )
descriptive_stats<-rbind(sumstats_epidemic_dummy, sumstats_epidemic_severity, sumstats_pop_dens_log, sumstats_night_lights_log, sumstats_conflict_intens_log, sumstats_political_exclud_group,sumstats_spatial_lag_social_unrest, sumstats_temporal_lag_social_unrest,sumstats_temporal_lag_repression, sumstats_democracy, sumstats_gdppc_log, sumstats_gdp_growth, sumstats_pop_log, sumstats_federalism, sumstats_hiv_log, sumstats_press_freedom)
descriptive_stats<-as.data.frame(descriptive_stats)
names(descriptive_stats) <- gsub("_epidemic.dummy", "", names(descriptive_stats))
descriptive_stats<-round(descriptive_stats, digits=3)
rownames(descriptive_stats) <- c("Epidemic Binary","Epidemic Severity","Population Density (log)","Stable Night Lights (log)","Conflict Intensity (log)", "Politically Excluded Group", "Spatially Weighted Lag Social Unrest" , "Social Unrest (t-1)", "Repression (t-1)","Democracy", "GDPPC (log)", "GDP Growth", "Population (log)", "Federalism", "HIV Cases (log)","Press Freedom")

# Print Descriptive Stats 
descriptive_stats

######################################################################
# Table A2 Imbalance Statistics, Pre-matched and Post-matched Sample
######################################################################

# Mean Diff (Pre)
nona<-pre_matched[,c("epidemic_dummy","pop_dens_log","night_lights_log","conflict_intens_log","democracy","gdppc_log","pop_log","physicians_per_1000")]
nona<-nona[complete.cases(nona), ]
pre_treat<-nona[nona$epidemic_dummy==1,]
pre_control<-nona[nona$epidemic_dummy==0,]
pre_treat<-colMeans(pre_treat[ names(pre_treat)])
pre_control<-colMeans(pre_control[ names(pre_control)])
pre_mean_diff<-as.data.frame(pre_control-pre_treat)[2:8,]

# L1 (Pre)
pre_imbalance <- imbalance(group=nona$epidemic_dummy,data=nona,drop=c("epidemic_dummy"))
pre_l1<-pre_imbalance$tab[3]

# Mean Diff (Post)
matched<-post_matched[,c("epidemic_dummy","pop_dens_log","night_lights_log","conflict_intens_log","democracy","gdppc_log","pop_log","physicians_per_1000")]
post_treat<-matched[matched$epidemic_dummy==1,]
post_control<-matched[matched$epidemic_dummy==0,]
post_treat<-colMeans(post_treat[ names(post_treat)])
post_control<-colMeans(post_control[ names(post_control)])
post_mean_diff<-as.data.frame(post_control-post_treat)[2:8,]

# L1 (Post)
post_imbalance <- imbalance(group=matched$epidemic_dummy,data=matched,drop=c("epidemic_dummy"))
post_l1<-post_imbalance$tab[3]

# Create table
imbalance_table<-round(cbind(pre_mean_diff, post_mean_diff, pre_l1, post_l1),digits=3)

# Print imbalance statistics table
imbalance_table

# "Overall L1 multivariate imbalance measure decreases from 0.96 for the pre-matched sample to 0.80 for the post-matched sample"
pre_overall_l1<-round(pre_imbalance$L1$L1, digits=2)
post_overall_l1<-round(post_imbalance$L1$L1, digits=2)

# Print
pre_overall_l1
post_overall_l1

#######################################################################################
# Table A3 - Multi-level Negative Binomial Regression Results, Epidemic Binary Variable
#######################################################################################

# Model 1 Pre-matched sample
taba3_m1 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=pre_matched)

# Model 2 Pre-matched sample
taba3_m2 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=pre_matched)

# Model 3 Post-matched sample
taba3_m3 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=post_matched)

# Model 4 Post-matched sample
taba3_m4 <- glmmTMB(social_unrest_count ~ epidemic_dummy + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=post_matched)

# Extract results
source(system.file("other_methods","extract.R",package="glmmTMB"))
extract.taba3_m1<-extract.glmmTMB(taba3_m1)
extract.taba3_m2<-extract.glmmTMB(taba3_m2)
extract.taba3_m3<-extract.glmmTMB(taba3_m3)
extract.taba3_m4<-extract.glmmTMB(taba3_m4)

# One-tailed test
extract.taba3_m1@pvalues<-(extract.taba3_m1@pvalues)/2
extract.taba3_m2@pvalues<-(extract.taba3_m2@pvalues)/2
extract.taba3_m3@pvalues<-(extract.taba3_m3@pvalues)/2
extract.taba3_m4@pvalues<-(extract.taba3_m4@pvalues)/2

# Print Regression Table
screenreg(list(extract.taba3_m1, extract.taba3_m2, extract.taba3_m3, extract.taba3_m4), custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"), stars = c(0.01, 0.05, 0.1))

#######################
# Model fit stats
#######################

# Observations
n_taba3_m1<-nrow(taba3_m1$frame)
n_taba3_m2<-nrow(taba3_m2$frame)
n_taba3_m3<-nrow(taba3_m3$frame)
n_taba3_m4<-nrow(taba3_m4$frame)

# Log Likelihood
log_taba3_m1<-logLik(taba3_m1)
log_taba3_m2<-logLik(taba3_m2)
log_taba3_m3<-logLik(taba3_m3)
log_taba3_m4<-logLik(taba3_m4)

# AIC
aic_taba3_m1<-AIC(taba3_m1)
aic_taba3_m2<-AIC(taba3_m2)
aic_taba3_m3<-AIC(taba3_m3)
aic_taba3_m4<-AIC(taba3_m4)

# BIC
bic_taba3_m1<-BIC(taba3_m1)
bic_taba3_m2<-BIC(taba3_m2)
bic_taba3_m3<-BIC(taba3_m3)
bic_taba3_m4<-BIC(taba3_m4)

# Create Model Fit Stats Table
all_taba3_m1<-rbind(n_taba3_m1,log_taba3_m1,aic_taba3_m1,bic_taba3_m1)
all_taba3_m2<-rbind(n_taba3_m2,log_taba3_m2,aic_taba3_m2,bic_taba3_m2)
all_taba3_m3<-rbind(n_taba3_m3,log_taba3_m3,aic_taba3_m3,bic_taba3_m3)
all_taba3_m4<-rbind(n_taba3_m4,log_taba3_m4,aic_taba3_m4,bic_taba3_m4)

# Print Model Fit Stats
all_taba3<-round(cbind(all_taba3_m1,all_taba3_m2,all_taba3_m3,all_taba3_m4),digits=2)
all_taba3

####################################################################################################

##########################################################################################
# Table A4 - Multi-level Negative Binomial Regression Results, Epidemic Severity Variable
##########################################################################################

# Model 1 Pre-matched sample
taba4_m1 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=pre_matched)

# Model 2 Pre-matched sample
taba4_m2 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=pre_matched)

# Model 3 Post-matched sample
taba4_m3 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=post_matched)

# Model 4 Post-matched sample
taba4_m4 <- glmmTMB(social_unrest_count ~ epidemic_severity + pop_dens_log + night_lights_log + conflict_intens_log + political_exclud_group + spatial_lag_social_unrest + temporal_lag_social_unrest + temporal_lag_repression + democracy + gdppc_log + gdp_growth + pop_log + federalism + hiv_log + press_freedom + (1 | GID_0) + (1 | GID_1), family = nbinom1, data=post_matched)

# Extract results
source(system.file("other_methods","extract.R",package="glmmTMB"))
extract.taba4_m1<-extract.glmmTMB(taba4_m1)
extract.taba4_m2<-extract.glmmTMB(taba4_m2)
extract.taba4_m3<-extract.glmmTMB(taba4_m3)
extract.taba4_m4<-extract.glmmTMB(taba4_m4)

# One-tailed test
extract.taba4_m1@pvalues<-(extract.taba4_m1@pvalues)/2
extract.taba4_m2@pvalues<-(extract.taba4_m2@pvalues)/2
extract.taba4_m3@pvalues<-(extract.taba4_m3@pvalues)/2
extract.taba4_m4@pvalues<-(extract.taba4_m4@pvalues)/2

# Print Regression Table
screenreg(list(extract.taba4_m1, extract.taba4_m2, extract.taba4_m3, extract.taba4_m4), custom.model.names = c("Model 1", "Model 2", "Model 3", "Model 4"), stars = c(0.01, 0.05, 0.1))

#######################
# Model fit stats
#######################

# Observations
n_taba4_m1<-nrow(taba4_m1$frame)
n_taba4_m2<-nrow(taba4_m2$frame)
n_taba4_m3<-nrow(taba4_m3$frame)
n_taba4_m4<-nrow(taba4_m4$frame)

# Log Likelihood
log_taba4_m1<-logLik(taba4_m1)
log_taba4_m2<-logLik(taba4_m2)
log_taba4_m3<-logLik(taba4_m3)
log_taba4_m4<-logLik(taba4_m4)

# AIC
aic_taba4_m1<-AIC(taba4_m1)
aic_taba4_m2<-AIC(taba4_m2)
aic_taba4_m3<-AIC(taba4_m3)
aic_taba4_m4<-AIC(taba4_m4)

# BIC
bic_taba4_m1<-BIC(taba4_m1)
bic_taba4_m2<-BIC(taba4_m2)
bic_taba4_m3<-BIC(taba4_m3)
bic_taba4_m4<-BIC(taba4_m4)

# Create Model Fit Stats Table
all_taba4_m1<-rbind(n_taba4_m1,log_taba4_m1,aic_taba4_m1,bic_taba4_m1)
all_taba4_m2<-rbind(n_taba4_m2,log_taba4_m2,aic_taba4_m2,bic_taba4_m2)
all_taba4_m3<-rbind(n_taba4_m3,log_taba4_m3,aic_taba4_m3,bic_taba4_m3)
all_taba4_m4<-rbind(n_taba4_m4,log_taba4_m4,aic_taba4_m4,bic_taba4_m4)

# Print Model Fit Stats
all_taba4<-round(cbind(all_taba4_m1,all_taba4_m2,all_taba4_m3,all_taba4_m4),digits=2)
all_taba4
